home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / qmath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  19.7 KB  |  1,185 lines  |  [TEXT/????]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision rational arithmetic primitive routines
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include "xmath.h"
  11.  
  12.  
  13. NUMBER _qzero_ =    { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  14. NUMBER _qone_ =        { { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  15. static NUMBER _qtwo_ =    { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  16. static NUMBER _qten_ =    { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
  17. NUMBER _qnegone_ =    { { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1 };
  18. NUMBER _qonehalf_ =    { { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1 };
  19.  
  20.  
  21. #if 0
  22. static char *abortmsg = "Calculation aborted";
  23. #endif
  24. static char *memmsg = "Not enough memory";
  25.  
  26.  
  27. /*
  28.  * Create another copy of a number.
  29.  *    q2 = qcopy(q1);
  30.  */
  31. NUMBER *
  32. qcopy(q)
  33.     register NUMBER *q;
  34. {
  35.     register NUMBER *r;
  36.  
  37.     r = qalloc();
  38.     r->num.sign = q->num.sign;
  39.     if (!isunit(q->num)) {
  40.         r->num.len = q->num.len;
  41.         r->num.v = alloc(r->num.len);
  42.         copyval(q->num, r->num);
  43.     }
  44.     if (!isunit(q->den)) {
  45.         r->den.len = q->den.len;
  46.         r->den.v = alloc(r->den.len);
  47.         copyval(q->den, r->den);
  48.     }
  49.     return r;
  50. }
  51.  
  52.  
  53. /*
  54.  * Convert a number to a normal integer.
  55.  *    i = qtoi(q);
  56.  */
  57. long
  58. qtoi(q)
  59.     register NUMBER *q;
  60. {
  61.     long i;
  62.     ZVALUE res;
  63.  
  64.     if (qisint(q))
  65.         return ztoi(q->num);
  66.     zquo(q->num, q->den, &res);
  67.     i = ztoi(res);
  68.     freeh(res.v);
  69.     return i;
  70. }
  71.  
  72.  
  73. /*
  74.  * Convert a normal integer into a number.
  75.  *    q = itoq(i);
  76.  */
  77. NUMBER *
  78. itoq(i)
  79.     long i;
  80. {
  81.     register NUMBER *q;
  82.  
  83.     if ((i >= -1) && (i <= 10)) {
  84.         switch ((int) i) {
  85.             case 0: q = &_qzero_; break;
  86.             case 1: q = &_qone_; break;
  87.             case 2: q = &_qtwo_; break;
  88.             case 10: q = &_qten_; break;
  89.             case -1: q = &_qnegone_; break;
  90.             default: q = NULL;
  91.         }
  92.         if (q)
  93.             return qlink(q);
  94.     }
  95.     q = qalloc();
  96.     itoz(i, &q->num);
  97.     return q;
  98. }
  99.  
  100.  
  101. /*
  102.  * Create a number from the given integral numerator and denominator.
  103.  *    q = iitoq(inum, iden);
  104.  */
  105. NUMBER *
  106. iitoq(inum, iden)
  107.     long inum, iden;
  108. {
  109.     register NUMBER *q;
  110.     long d;
  111.     BOOL sign;
  112.  
  113.     if (iden == 0)
  114.         error("Division by zero");
  115.     if (inum == 0)
  116.         return qlink(&_qzero_);
  117.     sign = 0;
  118.     if (inum < 0) {
  119.         sign = 1;
  120.         inum = -inum;
  121.     }
  122.     if (iden < 0) {
  123.         sign = 1 - sign;
  124.         iden = -iden;
  125.     }
  126.     d = iigcd(inum, iden);
  127.     inum /= d;
  128.     iden /= d;
  129.     if (iden == 1)
  130.         return itoq(sign ? -inum : inum);
  131.     q = qalloc();
  132.     if (inum != 1)
  133.         itoz(inum, &q->num);
  134.     itoz(iden, &q->den);
  135.     q->num.sign = sign;
  136.     return q;
  137. }
  138.  
  139.  
  140. /*
  141.  * Add two numbers to each other.
  142.  *    q3 = qadd(q1, q2);
  143.  */
  144. NUMBER *
  145. qadd(q1, q2)
  146.     register NUMBER *q1, *q2;
  147. {
  148.     NUMBER *r;
  149.     ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;
  150.  
  151.     r = qalloc();
  152.     /*
  153.      * If either number is an integer, then the result is easy.
  154.      */
  155.     if (qisint(q1) && qisint(q2)) {
  156.         zadd(q1->num, q2->num, &r->num);
  157.         return r;
  158.     }
  159.     if (qisint(q2)) {
  160.         zmul(q1->den, q2->num, &temp);
  161.         zadd(q1->num, temp, &r->num);
  162.         freeh(temp.v);
  163.         zcopy(q1->den, &r->den);
  164.         return r;
  165.     }
  166.     if (qisint(q1)) {
  167.         zmul(q2->den, q1->num, &temp);
  168.         zadd(q2->num, temp, &r->num);
  169.         freeh(temp.v);
  170.         zcopy(q2->den, &r->den);
  171.         return r;
  172.     }
  173.     /*
  174.      * Both arguments are true fractions, so we need more work.
  175.      * If the denominators are relatively prime, then the answer is the
  176.      * straightforward cross product result with no need for reduction.
  177.      */
  178.     zgcd(q1->den, q2->den, &d1);
  179.     if (isunit(d1)) {
  180.         freeh(d1.v);
  181.         zmul(q1->num, q2->den, &t1);
  182.         zmul(q1->den, q2->num, &t2);
  183.         zadd(t1, t2, &r->num);
  184.         freeh(t1.v);
  185.         freeh(t2.v);
  186.         zmul(q1->den, q2->den, &r->den);
  187.         return r;
  188.     }
  189.     /*
  190.      * The calculation is now more complicated.
  191.      * See Knuth Vol 2 for details.
  192.      */
  193.     zquo(q2->den, d1, &vpd1);
  194.     zquo(q1->den, d1, &upd1);
  195.     zmul(q1->num, vpd1, &t1);
  196.     zmul(q2->num, upd1, &t2);
  197.     zadd(t1, t2, &temp);
  198.     freeh(t1.v);
  199.     freeh(t2.v);
  200.     freeh(vpd1.v);
  201.     zgcd(temp, d1, &d2);
  202.     freeh(d1.v);
  203.     if (isunit(d2)) {
  204.         freeh(d2.v);
  205.         r->num = temp;
  206.         zmul(upd1, q2->den, &r->den);
  207.         freeh(upd1.v);
  208.         return r;
  209.     }
  210.     zquo(temp, d2, &r->num);
  211.     freeh(temp.v);
  212.     zquo(q2->den, d2, &temp);
  213.     freeh(d2.v);
  214.     zmul(temp, upd1, &r->den);
  215.     freeh(temp.v);
  216.     freeh(upd1.v);
  217.     return r;
  218. }
  219.  
  220.  
  221. /*
  222.  * Subtract one number from another.
  223.  *    q3 = qsub(q1, q2);
  224.  */
  225. NUMBER *
  226. qsub(q1, q2)
  227.     register NUMBER *q1, *q2;
  228. {
  229.     NUMBER t, *r;
  230.  
  231.     if (q1 == q2)
  232.         return qlink(&_qzero_);
  233.     if (qisint(q1) && qisint(q2)) {
  234.         r = qalloc();
  235.         zsub(q1->num, q2->num, &r->num);
  236.         return r;
  237.     }
  238.     t = *q2;
  239.     t.num.sign = !t.num.sign;
  240.     return qadd(q1, &t);
  241. }
  242.  
  243.  
  244. /*
  245.  * Increment a number by one.
  246.  */
  247. NUMBER *
  248. qinc(q)
  249.     NUMBER *q;
  250. {
  251.     NUMBER *r;
  252.  
  253.     r = qalloc();
  254.     if (qisint(q)) {
  255.         zadd(q->num, _one_, &r->num);
  256.         return r;
  257.     }
  258.     zadd(q->num, q->den, &r->num);
  259.     zcopy(q->den, &r->den);
  260.     return r;
  261. }
  262.  
  263.  
  264. /*
  265.  * Decrement a number by one.
  266.  */
  267. NUMBER *
  268. qdec(q)
  269.     NUMBER *q;
  270. {
  271.     NUMBER *r;
  272.  
  273.     r = qalloc();
  274.     if (qisint(q)) {
  275.         zsub(q->num, _one_, &r->num);
  276.         return r;
  277.     }
  278.     zsub(q->num, q->den, &r->num);
  279.     zcopy(q->den, &r->den);
  280.     return r;
  281. }
  282.  
  283.  
  284. #ifdef CODE
  285. /*
  286.  * Add a normal small integer value to an arbitrary number.
  287.  */
  288. NUMBER *
  289. qaddi(q1, n)
  290.     NUMBER *q1;
  291.     long n;
  292. {
  293.     NUMBER addnum;        /* temporary number */
  294.     HALF addval[2];        /* value of small number */
  295.     BOOL neg;        /* TRUE if number is neg */
  296.  
  297.     if (n == 0)
  298.         return qlink(q1);
  299.     if (n == 1)
  300.         return qinc(q1);
  301.     if (n == -1)
  302.         return qdec(q1);
  303.     if (qiszero(q1))
  304.         return itoq(n);
  305.     addnum.num.sign = 0;
  306.     addnum.num.len = 1;
  307.     addnum.num.v = addval;
  308.     addnum.den = _one_;
  309.     neg = (n < 0);
  310.     if (neg)
  311.         n = -n;
  312.     addval[0] = (HALF) n;
  313.     n = (((FULL) n) >> BASEB);
  314.     if (n) {
  315.         addval[1] = (HALF) n;
  316.         addnum.num.len = 2;
  317.     }
  318.     if (neg)
  319.         return qsub(q1, &addnum);
  320.     else
  321.         return qadd(q1, &addnum);
  322. }
  323. #endif
  324.  
  325.  
  326. /*
  327.  * Multiply two numbers.
  328.  *    q3 = qmul(q1, q2);
  329.  */
  330. NUMBER *
  331. qmul(q1, q2)
  332.     register NUMBER *q1, *q2;
  333. {
  334.     NUMBER *r;            /* returned value */
  335.     ZVALUE n1, n2, d1, d2;        /* numerators and denominators */
  336.     ZVALUE tmp;
  337.  
  338.     if (qiszero(q1) || qiszero(q2))
  339.         return qlink(&_qzero_);
  340.     if (qisone(q1))
  341.         return qlink(q2);
  342.     if (qisone(q2))
  343.         return qlink(q1);
  344.     if (qisint(q1) && qisint(q2)) {    /* easy results if integers */
  345.         r = qalloc();
  346.         zmul(q1->num, q2->num, &r->num);
  347.         return r;
  348.     }
  349.     n1 = q1->num;
  350.     n2 = q2->num;
  351.     d1 = q1->den;
  352.     d2 = q2->den;
  353.     if (iszero(d1) || iszero(d2))
  354.         error("Division by zero");
  355.     if (iszero(n1) || iszero(n2))
  356.         return qlink(&_qzero_);
  357.     if (!isunit(n1) && !isunit(d2)) {    /* possibly reduce */
  358.         zgcd(n1, d2, &tmp);
  359.         if (!isunit(tmp)) {
  360.             zquo(q1->num, tmp, &n1);
  361.             zquo(q2->den, tmp, &d2);
  362.         }
  363.         freeh(tmp.v);
  364.     }
  365.     if (!isunit(n2) && !isunit(d1)) {    /* again possibly reduce */
  366.         zgcd(n2, d1, &tmp);
  367.         if (!isunit(tmp)) {
  368.             zquo(q2->num, tmp, &n2);
  369.             zquo(q1->den, tmp, &d1);
  370.         }
  371.         freeh(tmp.v);
  372.     }
  373.     r = qalloc();
  374.     zmul(n1, n2, &r->num);
  375.     zmul(d1, d2, &r->den);
  376.     if (q1->num.v != n1.v)
  377.         freeh(n1.v);
  378.     if (q1->den.v != d1.v)
  379.         freeh(d1.v);
  380.     if (q2->num.v != n2.v)
  381.         freeh(n2.v);
  382.     if (q2->den.v != d2.v)
  383.         freeh(d2.v);
  384.     return r;
  385. }
  386.  
  387.  
  388. #if 0
  389. /*
  390.  * Multiply a number by a small integer.
  391.  *    q2 = qmuli(q1, n);
  392.  */
  393. NUMBER *
  394. qmuli(q, n)
  395.     NUMBER *q;
  396.     long n;
  397. {
  398.     NUMBER *r;
  399.     long d;            /* gcd of multiplier and denominator */
  400.     int sign;
  401.  
  402.     if ((n == 0) || qiszero(q))
  403.         return qlink(&_qzero_);
  404.     if (n == 1)
  405.         return qlink(q);
  406.     r = qalloc();
  407.     if (qisint(q)) {
  408.         zmuli(q->num, n, &r->num);
  409.         return r;
  410.     }
  411.     sign = 1;
  412.     if (n < 0) {
  413.         n = -n;
  414.         sign = -1;
  415.     }
  416.     d = zmodi(q->den, n);
  417.     d = iigcd(d, n);
  418.     zmuli(q->num, (n * sign) / d, &r->num);
  419.     (void) zdivi(q->den, d, &r->den);
  420.     return r;
  421. }
  422. #endif
  423.  
  424.  
  425. /*
  426.  * Divide two numbers (as fractions).
  427.  *    q3 = qdiv(q1, q2);
  428.  */
  429. NUMBER *
  430. qdiv(q1, q2)
  431.     register NUMBER *q1, *q2;
  432. {
  433.     NUMBER temp;
  434.  
  435.     if (q1 == q2) {
  436.         if (qiszero(q2))
  437.             error("Division by zero");
  438.         return qlink(&_qone_);
  439.     }
  440.     if (qisone(q1))
  441.         return qinv(q2);
  442.     temp.num = q2->den;
  443.     temp.den = q2->num;
  444.     temp.num.sign = temp.den.sign;
  445.     temp.den.sign = 0;
  446.     temp.links = 1;
  447.     return qmul(q1, &temp);
  448. }
  449.  
  450.  
  451. /*
  452.  * Divide a number by a small integer.
  453.  *    q2 = qdivi(q1, n);
  454.  */
  455. NUMBER *
  456. qdivi(q, n)
  457.     NUMBER *q;
  458.     long n;
  459. {
  460.     NUMBER *r;
  461.     long d;            /* gcd of divisor and numerator */
  462.     int sign;
  463.  
  464.     if (n == 0)
  465.         error("Division by zero");
  466.     if ((n == 1) || qiszero(q))
  467.         return qlink(q);
  468.     sign = 1;
  469.     if (n < 0) {
  470.         n = -n;
  471.         sign = -1;
  472.     }
  473.     r = qalloc();
  474.     d = zmodi(q->num, n);
  475.     d = iigcd(d, n);
  476.     (void) zdivi(q->num, d * sign, &r->num);
  477.     zmuli(q->den, n / d, &r->den);
  478.     return r;
  479. }
  480.  
  481.  
  482. /*
  483.  * Return the quotient when one number is divided by another.
  484.  * This works for fractional values also, and in all cases:
  485.  *    qquo(q1, q2) = int(q1 / q2).
  486.  */
  487. NUMBER *
  488. qquo(q1, q2)
  489.     register NUMBER *q1, *q2;
  490. {
  491.     ZVALUE num, den, res;
  492.     NUMBER *q;
  493.  
  494.     if (isunit(q1->num))
  495.         num = q2->den;
  496.     else if (isunit(q2->den))
  497.         num = q1->num;
  498.     else
  499.         zmul(q1->num, q2->den, &num);
  500.     if (isunit(q1->den))
  501.         den = q2->num;
  502.     else if (isunit(q2->num))
  503.         den = q1->den;
  504.     else
  505.         zmul(q1->den, q2->num, &den);
  506.     zquo(num, den, &res);
  507.     if ((num.v != q2->den.v) && (num.v != q1->num.v))
  508.         freeh(num.v);
  509.     if ((den.v != q2->num.v) && (den.v != q1->den.v))
  510.         freeh(den.v);
  511.     if (iszero(res)) {
  512.         freeh(res.v);
  513.         return qlink(&_qzero_);
  514.     }
  515.     if (isunit(res)) {
  516.         q = (res.sign ? &_qnegone_ : &_qone_);
  517.         freeh(res.v);
  518.         return qlink(q);
  519.     }
  520.     q = qalloc();
  521.     q->num = res;
  522.     return q;
  523. }
  524.  
  525.  
  526. /*
  527.  * Return the absolute value of a number.
  528.  *    q2 = qabs(q1);
  529.  */
  530. NUMBER *
  531. qabs(q)
  532.     register NUMBER *q;
  533. {
  534.     register NUMBER *r;
  535.  
  536.     if (q->num.sign == 0)
  537.         return qlink(q);
  538.     r = qalloc();
  539.     if (!isunit(q->num))
  540.         zcopy(q->num, &r->num);
  541.     if (!isunit(q->den))
  542.         zcopy(q->den, &r->den);
  543.     r->num.sign = 0;
  544.     return r;
  545. }
  546.  
  547.  
  548. /*
  549.  * Negate a number.
  550.  *    q2 = qneg(q1);
  551.  */
  552. NUMBER *
  553. qneg(q)
  554.     register NUMBER *q;
  555. {
  556.     register NUMBER *r;
  557.  
  558.     if (qiszero(q))
  559.         return qlink(&_qzero_);
  560.     r = qalloc();
  561.     if (!isunit(q->num))
  562.         zcopy(q->num, &r->num);
  563.     if (!isunit(q->den))
  564.         zcopy(q->den, &r->den);
  565.     r->num.sign = !q->num.sign;
  566.     return r;
  567. }
  568.  
  569.  
  570. /*
  571.  * Return the sign of a number (-1, 0 or 1)
  572.  */
  573. NUMBER *
  574. qsign(q)
  575.     NUMBER *q;
  576. {
  577.     if (qiszero(q))
  578.         return qlink(&_qzero_);
  579.     if (qisneg(q))
  580.         return qlink(&_qnegone_);
  581.     return qlink(&_qone_);
  582. }
  583.  
  584.  
  585. /*
  586.  * Invert a number.
  587.  *    q2 = qinv(q1);
  588.  */
  589. NUMBER *
  590. qinv(q)
  591.     register NUMBER *q;
  592. {
  593.     register NUMBER *r;
  594.  
  595.     if (qisunit(q)) {
  596.         r = (qisneg(q) ? &_qnegone_ : &_qone_);
  597.         return qlink(r);
  598.     }
  599.     if (qiszero(q))
  600.         error("Division by zero");
  601.     r = qalloc();
  602.     if (!isunit(q->num))
  603.         zcopy(q->num, &r->den);
  604.     if (!isunit(q->den))
  605.         zcopy(q->den, &r->num);
  606.     r->num.sign = q->num.sign;
  607.     r->den.sign = 0;
  608.     return r;
  609. }
  610.  
  611.  
  612. /*
  613.  * Return just the numerator of a number.
  614.  *    q2 = qnum(q1);
  615.  */
  616. NUMBER *
  617. qnum(q)
  618.     register NUMBER *q;
  619. {
  620.     register NUMBER *r;
  621.  
  622.     if (qisint(q))
  623.         return qlink(q);
  624.     if (isunit(q->num)) {
  625.         r = (qisneg(q) ? &_qnegone_ : &_qone_);
  626.         return qlink(r);
  627.     }
  628.     r = qalloc();
  629.     zcopy(q->num, &r->num);
  630.     return r;
  631. }
  632.  
  633.  
  634. /*
  635.  * Return just the denominator of a number.
  636.  *    q2 = qden(q1);
  637.  */
  638. NUMBER *
  639. qden(q)
  640.     register NUMBER *q;
  641. {
  642.     register NUMBER *r;
  643.  
  644.     if (qisint(q))
  645.         return qlink(&_qone_);
  646.     r = qalloc();
  647.     zcopy(q->den, &r->num);
  648.     return r;
  649. }
  650.  
  651.  
  652. /*
  653.  * Return the fractional part of a number.
  654.  *    q2 = qfrac(q1);
  655.  */
  656. NUMBER *
  657. qfrac(q)
  658.     register NUMBER *q;
  659. {
  660.     register NUMBER *r;
  661.  
  662.     if (qisint(q))
  663.         return qlink(&_qzero_);
  664.     if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
  665.         (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
  666.             return qlink(q);
  667.     r = qalloc();
  668.     zmod(q->num, q->den, &r->num);
  669.     zcopy(q->den, &r->den);
  670.     r->num.sign = q->num.sign;
  671.     return r;
  672. }
  673.  
  674.  
  675. /*
  676.  * Return the integral part of a number.
  677.  *    q2 = qint(q1);
  678.  */
  679. NUMBER *
  680. qint(q)
  681.     register NUMBER *q;
  682. {
  683.     register NUMBER *r;
  684.  
  685.     if (qisint(q))
  686.         return qlink(q);
  687.     if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
  688.         (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
  689.             return qlink(&_qzero_);
  690.     r = qalloc();
  691.     zquo(q->num, q->den, &r->num);
  692.     return r;
  693. }
  694.  
  695.  
  696. /*
  697.  * Compute the square of a number.
  698.  */
  699. NUMBER *
  700. qsquare(q)
  701.     register NUMBER *q;
  702. {
  703.     ZVALUE num, den;
  704.  
  705.     if (qiszero(q))
  706.         return qlink(&_qzero_);
  707.     if (qisunit(q))
  708.         return qlink(&_qone_);
  709.     num = q->num;
  710.     den = q->den;
  711.     q = qalloc();
  712.     if (!isunit(num))
  713.         zsquare(num, &q->num);
  714.     if (!isunit(den))
  715.         zsquare(den, &q->den);
  716.     return q;
  717. }
  718.  
  719.  
  720. /*
  721.  * Shift an integer by a given number of bits. This multiplies the number
  722.  * by the appropriate power of two.  Positive numbers shift left, negative
  723.  * ones shift right.  Low bits are truncated when shifting right.
  724.  */
  725. NUMBER *
  726. qshift(q, n)
  727.     NUMBER *q;
  728.     long n;
  729. {
  730.     register NUMBER *r;
  731.  
  732.     if (qisfrac(q))
  733.         error("Shift of non-integer");
  734.     if (qiszero(q) || (n == 0))
  735.         return qlink(q);
  736.     if (n <= -(q->num.len * BASEB))
  737.         return qlink(&_qzero_);
  738.     r = qalloc();
  739.     zshift(q->num, n, &r->num);
  740.     return r;
  741. }
  742.  
  743.  
  744. /*
  745.  * Scale a number by a power of two, as in:
  746.  *    ans = q * 2^n.
  747.  * This is similar to shifting, except that fractions work.
  748.  */
  749. NUMBER *
  750. qscale(q, pow)
  751.     NUMBER *q;
  752.     long pow;
  753. {
  754.     long numshift, denshift, tmp;
  755.     NUMBER *r;
  756.  
  757.     if (qiszero(q) || (pow == 0))
  758.         return qlink(q);
  759.     if ((pow > 1000000L) || (pow < -1000000L))
  760.         error("Very large scale value");
  761.     numshift = isodd(q->num) ? 0 : zlowbit(q->num);
  762.     denshift = isodd(q->den) ? 0 : zlowbit(q->den);
  763.     if (pow > 0) {
  764.         tmp = pow;
  765.         if (tmp > denshift)
  766.         tmp = denshift;
  767.         denshift = -tmp;
  768.         numshift = (pow - tmp);
  769.     } else {
  770.         pow = -pow;
  771.         tmp = pow;
  772.         if (tmp > numshift)
  773.             tmp = numshift;
  774.         numshift = -tmp;
  775.         denshift = (pow - tmp);
  776.     }
  777.     r = qalloc();
  778.     if (numshift)
  779.         zshift(q->num, numshift, &r->num);
  780.     else
  781.         zcopy(q->num, &r->num);
  782.     if (denshift)
  783.         zshift(q->den, denshift, &r->den);
  784.     else
  785.         zcopy(q->den, &r->den);
  786.     return r;
  787. }
  788.  
  789.  
  790. /*
  791.  * Return the minimum of two numbers.
  792.  */
  793. NUMBER *
  794. qmin(q1, q2)
  795.     NUMBER *q1, *q2;
  796. {
  797.     if (qrel(q1, q2) > 0)
  798.         q1 = q2;
  799.     return qlink(q1);
  800. }
  801.  
  802.  
  803. /*
  804.  * Return the maximum of two numbers.
  805.  */
  806. NUMBER *
  807. qmax(q1, q2)
  808.     NUMBER *q1, *q2;
  809. {
  810.     if (qrel(q1, q2) < 0)
  811.         q1 = q2;
  812.     return qlink(q1);
  813. }
  814.  
  815.  
  816. /*
  817.  * Perform the logical OR of two integers.
  818.  */
  819. NUMBER *
  820. qor(q1, q2)
  821.     NUMBER *q1, *q2;
  822. {
  823.     register NUMBER *r;
  824.  
  825.     if (qisfrac(q1) || qisfrac(q2))
  826.         error("Non-integers for logical or");
  827.     if ((q1 == q2) || qiszero(q2))
  828.         return qlink(q1);
  829.     if (qiszero(q1))
  830.         return qlink(q2);
  831.     r = qalloc();
  832.     zor(q1->num, q2->num, &r->num);
  833.     return r;
  834. }
  835.  
  836.  
  837. /*
  838.  * Perform the logical AND of two integers.
  839.  */
  840. NUMBER *
  841. qand(q1, q2)
  842.     NUMBER *q1, *q2;
  843. {
  844.     register NUMBER *r;
  845.     ZVALUE res;
  846.  
  847.     if (qisfrac(q1) || qisfrac(q2))
  848.         error("Non-integers for logical and");
  849.     if (q1 == q2)
  850.         return qlink(q1);
  851.     if (qiszero(q1) || qiszero(q2))
  852.         return qlink(&_qzero_);
  853.     zand(q1->num, q2->num, &res);
  854.     if (iszero(res)) {
  855.         freeh(res.v);
  856.         return qlink(&_qzero_);
  857.     }
  858.     r = qalloc();
  859.     r->num = res;
  860.     return r;
  861. }
  862.  
  863.  
  864. /*
  865.  * Perform the logical XOR of two integers.
  866.  */
  867. NUMBER *
  868. qxor(q1, q2)
  869.     NUMBER *q1, *q2;
  870. {
  871.     register NUMBER *r;
  872.     ZVALUE res;
  873.  
  874.     if (qisfrac(q1) || qisfrac(q2))
  875.         error("Non-integers for logical xor");
  876.     if (q1 == q2)
  877.         return qlink(&_qzero_);
  878.     if (qiszero(q1))
  879.         return qlink(q2);
  880.     if (qiszero(q2))
  881.         return qlink(q1);
  882.     zxor(q1->num, q2->num, &res);
  883.     if (iszero(res)) {
  884.         freeh(res.v);
  885.         return qlink(&_qzero_);
  886.     }
  887.     r = qalloc();
  888.     r->num = res;
  889.     return r;
  890. }
  891.  
  892.  
  893. #if 0
  894. /*
  895.  * Return the number whose binary representation only has the specified
  896.  * bit set (counting from zero).  This thus produces a given power of two.
  897.  */
  898. NUMBER *
  899. qbitvalue(n)
  900.     long n;
  901. {
  902.     register NUMBER *r;
  903.  
  904.     if (n <= 0)
  905.         return qlink(&_qone_);
  906.     r = qalloc();
  907.     zbitvalue(n, &r->num);
  908.     return r;
  909. }
  910.  
  911.  
  912. /*
  913.  * Test to see if the specified bit of a number is on (counted from zero).
  914.  * Returns TRUE if the bit is set, or FALSE if it is not.
  915.  *    i = qbittest(q, n);
  916.  */
  917. BOOL
  918. qbittest(q, n)
  919.     register NUMBER *q;
  920.     long n;
  921. {
  922.     int x, y;
  923.  
  924.     if ((n < 0) || (n >= (q->num.len * BASEB)))
  925.         return FALSE;
  926.     x = q->num.v[n / BASEB];
  927.     y = (1 << (n % BASEB));
  928.     return ((x & y) != 0);
  929. }
  930. #endif
  931.  
  932.  
  933. /*
  934.  * Return the precision of a number (usually for examining an epsilon value).
  935.  * This is the largest power of two whose reciprocal is not smaller in absolute
  936.  * value than the specified number.  For example, qbitprec(1/100) = 6.
  937.  * Numbers larger than one have a precision of zero.
  938.  */
  939. long
  940. qprecision(q)
  941.     NUMBER *q;
  942. {
  943.     long r;
  944.  
  945.     if (qisint(q))
  946.         return 0;
  947.     if (isunit(q->num))
  948.         return zhighbit(q->den);
  949.     r = zhighbit(q->den) - zhighbit(q->num) - 1;
  950.     if (r < 0)
  951.         r = 0;
  952.     return r;
  953. }
  954.  
  955.  
  956. #if 0
  957. /*
  958.  * Return an integer indicating the sign of a number (-1, 0, or 1).
  959.  *    i = qtst(q);
  960.  */
  961. FLAG
  962. qtest(q)
  963.     register NUMBER *q;
  964. {
  965.     if (!ztest(q->num))
  966.         return 0;
  967.     if (q->num.sign)
  968.         return -1;
  969.     return 1;
  970. }
  971. #endif
  972.  
  973.  
  974. /*
  975.  * Compare two numbers and return an integer indicating their relative size.
  976.  *    i = qrel(q1, q2);
  977.  */
  978. FLAG
  979. qrel(q1, q2)
  980.     register NUMBER *q1, *q2;
  981. {
  982.     ZVALUE z1, z2;
  983.     long wc1, wc2;
  984.     int sign;
  985.     int z1f = 0, z2f = 0;
  986.  
  987.     if (q1 == q2)
  988.         return 0;
  989.     sign = q2->num.sign - q1->num.sign;
  990.     if (sign)
  991.         return sign;
  992.     if (qiszero(q2))
  993.         return !qiszero(q1);
  994.     if (qiszero(q1))
  995.         return -1;
  996.     /*
  997.      * Make a quick comparison by calculating the number of words resulting as
  998.      * if we multiplied through by the denominators, and then comparing the
  999.      * word counts.
  1000.      */
  1001.     sign = 1;
  1002.     if (qisneg(q1))
  1003.         sign = -1;
  1004.     wc1 = q1->num.len + q2->den.len;
  1005.     wc2 = q2->num.len + q1->den.len;
  1006.     if (wc1 < wc2 - 1)
  1007.         return -sign;
  1008.     if (wc2 < wc1 - 1)
  1009.         return sign;
  1010.     /*
  1011.      * Quick check failed, must actually do the full comparison.
  1012.      */
  1013.     if (isunit(q2->den))
  1014.         z1 = q1->num;
  1015.     else if (isone(q1->num))
  1016.         z1 = q2->den;
  1017.     else {
  1018.         z1f = 1;
  1019.         zmul(q1->num, q2->den, &z1);
  1020.     }
  1021.     if (isunit(q1->den))
  1022.         z2 = q2->num;
  1023.     else if (isone(q2->num))
  1024.         z2 = q1->den;
  1025.     else {
  1026.         z2f = 1;
  1027.         zmul(q2->num, q1->den, &z2);
  1028.     }
  1029.     sign = zrel(z1, z2);
  1030.     if (z1f)
  1031.         freeh(z1.v);
  1032.     if (z2f)
  1033.         freeh(z2.v);
  1034.     return sign;
  1035. }
  1036.  
  1037.  
  1038. /*
  1039.  * Compare two numbers to see if they are equal.
  1040.  * This differs from qrel in that the numbers are not ordered.
  1041.  * Returns TRUE if they differ.
  1042.  */
  1043. BOOL
  1044. qcmp(q1, q2)
  1045.     register NUMBER *q1, *q2;
  1046. {
  1047.     if (q1 == q2)
  1048.         return FALSE;
  1049.     if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) ||
  1050.         (q2->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) ||
  1051.         (*q1->den.v != *q2->den.v))
  1052.             return TRUE;
  1053.     if (zcmp(q1->num, q2->num))
  1054.         return TRUE;
  1055.     if (qisint(q1))
  1056.         return FALSE;
  1057.     return zcmp(q1->den, q2->den);
  1058. }
  1059.  
  1060.  
  1061. /*
  1062.  * Compare a number against a normal small integer.
  1063.  * Returns 1, 0, or -1, according to whether the first number is greater,
  1064.  * equal, or less than the second number.
  1065.  *    n = qreli(q, n);
  1066.  */
  1067. FLAG
  1068. qreli(q, n)
  1069.     NUMBER *q;
  1070.     long n;
  1071. {
  1072.     int sign;
  1073.     ZVALUE num;
  1074.     HALF h2[2];
  1075.     NUMBER q2;
  1076.  
  1077.     sign = ztest(q->num);        /* do trivial sign checks */
  1078.     if (sign == 0) {
  1079.         if (n > 0)
  1080.             return -1;
  1081.         return (n < 0);
  1082.     }
  1083.     if ((sign < 0) && (n >= 0))
  1084.         return -1;
  1085.     if ((sign > 0) && (n <= 0))
  1086.         return 1;
  1087.     n *= sign;
  1088.     if (n == 1) {            /* quick check against 1 or -1 */
  1089.         num = q->num;
  1090.         num.sign = 0;
  1091.         return (sign * zrel(num, q->den));
  1092.     }
  1093.     num.sign = (sign < 0);
  1094.     num.len = 1 + (n >= BASE);
  1095.     num.v = h2;
  1096.     h2[0] = (n & BASE1);
  1097.     h2[1] = (n >> BASEB);
  1098.     if (isunit(q->den))    /* integer compare if no denominator */
  1099.         return zrel(q->num, num);
  1100.     q2.num = num;
  1101.     q2.den = _one_;
  1102.     q2.links = 1;
  1103.     return qrel(q, &q2);    /* full fractional compare */
  1104. }
  1105.  
  1106.  
  1107. /*
  1108.  * Compare a number against a small integer to see if they are equal.
  1109.  * Returns TRUE if they differ.
  1110.  */
  1111. BOOL
  1112. qcmpi(q, n)
  1113.     NUMBER *q;
  1114.     long n;
  1115. {
  1116.     long len;
  1117.  
  1118.     len = q->num.len;
  1119.     if ((len > 2) || qisfrac(q) || (q->num.sign != (n < 0)))
  1120.         return TRUE;
  1121.     if (n < 0)
  1122.         n = -n;
  1123.     if (((HALF)(n)) != q->num.v[0])
  1124.         return TRUE;
  1125.     n = ((FULL) n) >> BASEB;
  1126.     return (((n != 0) != (len == 2)) || (n != q->num.v[1]));
  1127. }
  1128.  
  1129.  
  1130. /*
  1131.  * Number node allocation routines
  1132.  */
  1133.  
  1134. #define    NNALLOC    1000
  1135.  
  1136. union allocNode {
  1137.     NUMBER    num;
  1138.     union allocNode    *link;
  1139. };
  1140.  
  1141. static union allocNode    *freeNum;
  1142.  
  1143.  
  1144. NUMBER *
  1145. qalloc()
  1146. {
  1147.     register union allocNode *temp;
  1148.  
  1149.     if (!freeNum) {
  1150.         freeNum = (union allocNode *)
  1151.         malloc(sizeof (NUMBER) * NNALLOC);
  1152.         if (freeNum == NULL)
  1153.             error(memmsg);
  1154.         temp = freeNum;
  1155.         while (temp != freeNum + NNALLOC - 2) {
  1156.             temp->link = temp+1;
  1157.             ++temp;
  1158.         }
  1159.     }
  1160.     temp = freeNum;
  1161.     freeNum = temp->link;
  1162.     temp->num.links = 1;
  1163.     temp->num.num = _one_;
  1164.     temp->num.den = _one_;
  1165.     return &temp->num;
  1166. }
  1167.  
  1168.  
  1169. void
  1170. qfreenum(q)
  1171.     register NUMBER *q;
  1172. {
  1173.     union allocNode *a;
  1174.  
  1175.     if (q == NULL)
  1176.         return;
  1177.     freeh(q->num.v);
  1178.     freeh(q->den.v);
  1179.     a = (union allocNode *) q;
  1180.     a->link = freeNum;
  1181.     freeNum = a;
  1182. }
  1183.  
  1184. /* END CODE */
  1185.